Distinct mtDNA lineages in free‐ranging Ammotragus(aoudad) from the United States indicate multiple introductions from northern Africa

Abstract Translocation records indicate aoudad (Ammotragus lervia) populations in the United States are a product of multiple human‐mediated introductions. Two mitochondrial markers (cytochrome b, cytb; displacement loop, D loop) and one nuclear gene (prion protein gene exon 3, PRNP) were used to determine: (1) genetic variation, (2) if genetic units correspond to taxonomic designations, (3) the number and geographic origin of translocations, and (4) divergence times. Three phylogenetic analyses (Bayesian inference, maximum likelihood, and parsimony) produced similar topologies with two clades (I and II). Clade I contained progeny of individuals resulting from introductions to Texas and Spain, and individuals from Algeria. Individuals in Clade II were progeny of past introductions to the United States and Europe, and northern Algeria. Clade II was subdivided into two subclades (A and B) representing two haplogroups. No genetic variation was detected in the PRNP sequences. Three haplogroups appeared to correspond to the subspecies A. l. lervia and A. l. sahariensis whose native distribution includes northwestern Africa. Network analyses assigned haplogroups to two major groups similar to those depicted in the phylogenetic analyses. Genetic distances ranged from 0.80% to 5.17% and 2.99% to 15.42% for cytb and D loop, respectively; and were higher than normally recovered for caprids, warranting a reexamination of subspecific status. Divergence dates indicated a major split between A. l. lervia and A. l. sahariensis circa 2.38 mya. Together, the high level of genetic divergences among US populations and apparent presence of two subspecies of aoudad in the United States support the hypothesis of multiple introductions from multiple sources.

l. blainei (Kordofan Aoudad), A. l. fassini (Libyan Aoudad), A. l. angusi (Aïr Aoudad), and A. l. sahariensis (Saharan Aoudad). Little is known concerning genetic variation and phylogeographic patterns of diversity among these putative subspecies. In reviewing the taxonomy of Caprini, Groves and Grubb (2011) questioned whether A. lervia was a single species. The sole genetic study (Derouiche et al., 2020) included wild-caught individuals from the Algerian provinces of Béchar, Illizi, and Tamanrasset; as well as semi-captive individuals obtained from several preserves and zoos in Algeria. Their results based on mitochondrial DNA sequences, reflected a Mediterranean and Saharan divergence which seems to correspond to the two subspecies (A. l. sahariensis and A. l. lervia). To date, no genetic information is available for the other four subspecies, complicating the resolution of phylogenetic relationships among subspecies within their native range, as well as their role in determining source-stock origins of introduced populations throughout the world.
A recent study (Stipoljev et al., 2021) used a mitochondrial marker (displacement loop, D loop) and microsatellite loci to determine genetic diversity and population structure in introduced populations of aoudad in Croatia, Czech Republic, and Spain. They reported evidence of four haplotypes and based on nuclear data they identified significant structure among populations. Stipoljev et al. (2021) suggested that the Almería haplotype probably was associated with A. l. sahariensis, whereas the other three haplotypes were of admixed origin presumably assignable to A. l. lervia. Stipoljev et al. (2021) reported low genetic diversity among populations (low number of detectable alleles and high number of shared alleles), consistent with a history of recent introductions (<50 years) and a small number of founding individuals from documented translocations.
Although aoudad are listed as "vulnerable" in their native range by the IUCN Red List of Threatened Species (Cassinello et al., 2008), substantial populations have been established in Europe and the southwestern United States (California, New Mexico, and Texas).
In fact, the number of non-native aoudad in the United States are thought to outnumber those existing in the native range (Cassinello et al., 2008;Stipoljev et al., 2021). Based on zoo records, it appears that aoudad initially were imported into the New York Zoological Park and the National Zoological Park in the United States, circa 1900 (Mungall & Sheffield, 1994;Ogren, 1959). Mexico (Barrett, 1980;Mungall & Sheffield, 1994;Ogren, 1965). In  (DeArment, 1971, Mungall & Sheffield, 1994Ogren, 1965). Further, throughout the 1950s and 1970s, private ranches independently introduced aoudad (Simpson & Krysl, 1981) into the eastern, central, and southwestern portions of Texas (Mungall & Sheffield, 1994). It is unclear whether these translocations included individuals previously established in Texas or were products of additional importations from their native range, zoos, or introduced populations in Europe.
At present, similar translocation efforts and population expansion continue with aoudad now being common throughout the western two-thirds of Texas. Currently, >30,000 free-ranging aoudad are estimated to occur in Texas; with most populations residing in the Trans-Pecos region, followed by the Edwards Plateau and Panhandle regions (F. Hernández, TPWD, personal communication; Traweek & Welch 1992; Figure 3); although aoudad occur in other ecoregions F I G U R E 1 Photograph of Ammotragus lervia in Palo Duro Canyon, Texas as a result of private introductions in high-fenced, non-free-ranging operations, and subsequent escapees (Schmidly & Bradley, 2016).
Beyond the aforementioned zoo records and discussions pre-

sented in the Symposium on Ecology and Management of Barbary
Sheep (Simpson, 1980) and New Mexico Game and Fish Bulletins (Ogren, 1962(Ogren, , 1965, there is little information relative to sourcestock origins and introductions into the United States and phylogenetic relationships of the six-native subspecies in northern Africa. Most experts (Gray, 1985;Ogren, 1965) surmised that A. l. lervia served as the source-stock for North American zoo introductions.
Although most zoos did not retain detailed source of origin records (geographic history) in the early 1900s, those that documented the movement/transfers of aoudad between zoos (e.g., trading between the National Zoological Park and New York Zoological Park) were not always in agreement. Further, given that private ranches were not obligated to maintain source-stock information or geographic history records, they generally did not provide point of origin data for imported aoudad. Consequently, the multiple imports and subsequent purchase or trading of zoo progeny by private individuals, coupled with a paucity of source/origin documentation, make it difficult to discern if more than one subspecies of aoudad was introduced into the United States during this time frame.
Herein, we aim to provide the first broad-scale geographic examination of aoudad in the southwestern United States (California, New Mexico, and Texas). The goals of this study were to: (1) assess genetic variation in free-ranging populations of aoudad in Texas, California, and New Mexico, (2) determine if genetic units correspond to existing taxonomic designations, (3) ascertain the number and geographic source of introductions, and (4) provide estimates of approximate divergence times among haplogroups. Mitochondrial markers, cytochrome b (cytb), and displacement loop (D loop) and one nuclear gene (prion protein gene exon 3, PRNP) were used to determine the magnitude of genetic variation. The cytb marker is used widely as a proxy to measure genetic divergence among species and subspecies (Baker & Bradley, 2006;Bradley & Baker, 2001;Larsen et al., 2010). The D loop marker was selected because its rapid rate of nucleotide sequence evolution makes it ideal for examining genetic variation between and within populations and as a measure of nucleotide and haplotype diversity and other genetic indices (Latch et al., 2009;Mendez-Harclerode et al., 2007;Stipoljev et al., 2021).
DNA sequences for PRNP were available from other ongoing studies on aoudad and bighorn sheep; consequently, it was included as a means for detecting genetic variation in the nuclear genome.

| Sampling
A total of 232 aoudad samples were collected between 2018 and 2021, these included 209 samples from free-ranging individuals, 19 samples from pedigreed captive individuals from the Fossil Rim F I G U R E 2 Map depicting the distribution of aoudad in North Africa based on Cassinello et al. (2008) and Cassinello (2013). Populations from northeastern Chad have been assigned to A. l. blainei (Alados et al., 1988) and A. l. sahariensis (Cassinello, 2013), in which resolving this is beyond the scope of this study. Therefore, this population is indicated by cross hashing to reflect its uncertainty Wildlife Center (Texas), and four historic samples (three from Palo Pinto County, Texas circa 1985 and one from Garza County, Texas circa 1991), and were used in this study (see Figure 3; Appendix).
Samples (ear clip, muscle, liver, whole blood, and/or dried muscle) obtained in the United States were acquired through five meth- Collection. Tissue samples obtained were either: (1) stored on ice and eventually frozen at −20°C or (2) immediately flash-frozen in liquid nitrogen. All tissue samples and museum specimens were deposited into the NSRL. Specimens collected in the above procedures followed methods outlined in the guidelines of the American Society of Mammalogists (Sikes et al., 2016) and protocols approved by the Texas Tech University Animal Care and Use Committee (protocols #17023-02 and 20002-01).

| DNA sequencing
Genomic DNA (gDNA) was extracted from 0.1 g ear clip, muscle, liver, or 100 µl blood (stored in standard collection tubes containing EDTA) using the Qiagen DNeasy blood and tissue extraction kit (Qiagen, Valencia, California). The entire cytb gene (1,143 bp)

A. l. lervia (Clade II-A)
A. l. lervia (Clade II-B) was amplified using the polymerase chain reaction (PCR) method (Saiki et al., 1988) with primers LGL765 (forward, Bickham et al., 1995) and LGL766 (reverse, Bickham et al., 2004) LGL765 (Bickham et al., 1995(Bickham et al., , 2004, 870R (Peppers et al., 2002), and F1 (Whiting et al., 2003). Cycle sequencing products subsequently were purified using Sephadex filtration (GE Healthcare) and centrifugation methods, followed by dehydration. Purified sequencing products were analyzed on an ABI 3730xl automated sequencer  Data obtained from the cytb and D loop datasets were used to direct selection of individuals to be examined for the prion protein exon 3 gene (PRNP). Primers used to amplify the complete PRNP gene (771 bp) were MD582F (forward, Jewell et al., 2005) and MD1479RC (reverse, Jewell et al., 2005). Thermal profiles for PCR were as follows: a hot start of 80°C, initial denaturation at 95°C for 2 min, followed by 35 cycles of denaturation at 95°C for 30 s, annealing at 54°C for 45 s, and extension at 72°C for 1 min, with a final extension at 72°C for 15 min. Primers used to cycle sequence the products included MD582F, MD1479RC, 12, and 3FL1 (Jewell et al., 2005). All cytb, D loop, and PRNP sequences obtained in this study were deposited in GenBank (accession numbers: Additional sequence data for both mitochondrial markers (cytb:

| Phylogenetic analyses
In the following analyses, data were obtained from three independent studies: cytb only (Derouiche et al., 2020), D loop only (Stipoljev et al., 2021), and cytb and D loop combined (this study A parsimony analysis (PAUP* version 4.0a169, Swofford, 2003) was conducted to identify synapomorphies indicative of taxonomic identifications. Parsimony characters were assigned equal weight and variable nucleotide positions were treated as unordered, discrete characters with four possible states: A, C, G, and T. Phylogenetically uninformative characters were removed from the analysis. The most-parsimonious trees were estimated using the heuristic search and tree-bisection-reconnection option. A strict consensus tree was generated from the population of most-parsimonious trees and a subsequent bootstrap analysis (Felsenstein, 1985) with 1,000 iterations and the "fast" step-wise option selected to evaluate nodal support.
Eighty-eight maximum likelihood (ML) models were evaluated

| Genetic divergence
Genetic distance values for selected taxa and mitochondrial haplogroups were estimated using the Kimura 2-parameter model of evolution (Kimura, 1980) and the Tamura-Nei model of evolution (Tamura & Nei, 1993) for the cytb and D loop datasets, respectively, using the program MEGA X (Kumar et al., 2018). The resulting values calculated from the mitochondrial markers were used to examine levels of genetic divergence pertaining to the genetic species concept outlined in Bradley and Baker (2001) and Baker and Bradley (2006). Sequences of both cytb and D loop for closely related taxa (Yang et al., 2013) were obtained from GenBank to provide comparative genetic distance values.

| Divergence dating
A Molecular Clock Test (ML, Kumar et al., 2018) was used to accept or reject a strict molecular clock. This result was used with the software program BEAST v2.6.2 (Bouckaert et al., 2019) as a prior to estimate molecular timelines associated with phylogenetic divergence.
Divergence dates for aoudad were estimated from the reduced cytb dataset (as explained above) using Rupicapra rupicapra, R. pyrenaica, and Arabitragus jayakari as outgroup taxa. Fossil calibrations were x 10 7 generations was analyzed with log and tree files, which were then combined to generate divergence date estimates and produce a maximum clade credibility tree. The program Tracer (Rambaut et al., 2018) was used to examine number of successful MCMC iterations, stability of topological convergence, and Effective Sample Size using a value >200 as indicative of a minimal threshold for all parameters.
The program TreeAnnotator (Bouckaert et al., 2019) subsequently was used to obtain an estimate of the final phylogenetic tree with divergence dates assigned to nodes.

| Characterization of PRNP
PRNPsequences were proofed using Sequencher 4.10.1 software (Gene Codes Corporation, Ann Arbor, Michigan) and heterozygous nucleotide base positions were visually determined using chromatograms. The program MEGA X (Kumar et al., 2018) was used to translate the nucleotide sequences to protein, allowing for the detection of any non-synonymous substitutions based on an outgroup comparison to closely related genera (Capra and Ovis). In particular, three codons known to be of importance (Goldmann, 2008) in identifying susceptibility to scrapie in domestic sheep and goats were examined. These polymorphic codons included: A136V/T, R154H/L, and Q171R/H/K (referenced using traditional amino acid terminology, Dunnen & Antonarakis, 2000).

| Cytochrome b dataset
For the cytb dataset, the preliminary neighbor-joining analysis (not shown) of 249 individuals was used to assign haplotype affiliation of all individuals to the sampled localities (see Figure 3). From this analysis, a reduced dataset (n = 23) was obtained when identical sequences were removed. This final dataset was used for all subsequent analyses involving phylogenetics, genetic distances, and genetic indices. The three phylogenetic analyses (BI, ML, and parsimony) generated similar topologies in the cytb dataset; consequently, each analysis is discussed in detail below; however, only the topology obtained from the BI analysis is shown (Figure 4). Although there was substantial variation among individuals in terminal nodes, these associations were collapsed due to lack of nodal support.
In the BI analysis, two supported clades were identified (I and II

| Insertion and deletion events
Ten insertion and deletions events (indels), representing a total of 64 nucleotide substitutions, were detected in the D loop dataset (Table 1). These indels occurred between nucleotides sites 15,556 to 15,772 (aoudad reference genome, NC009510). Deletions ranged from a single nucleotide to 16 bp, whereas insertions ranged from two to 15 nucleotides. Several of the indels were of phylogenetic relevance, for example, the first deletion event was 15 bp and restricted to individuals in Clade II, whereas the second deletion event was 16 bp in length and restricted to individuals in Clade I.
Collectively, these indels contributed to the greater branch lengths depicted in the D loop topology relative to the cytb dataset.

| Genetic distances
Estimation of Kimura-2 parameter (Kimura, 1980) genetic distances ( II-A; and 5.12% % between Clade I and Subclade II-B (Figure 4). Table 3) obtained from the D loop dataset were estimated using the Tamura and Nei model of evolution (Tamura & Nei, 1993 II-B ( Figure 5). TA B L E 2 Average genetic distances of cytb sequences estimated using the Kimura 2-parameter model of evolution (Kimura, 1980) for selected comparisons of aoudad and taxa of Family Bovidae

The major split between individuals assigned to Subclade II-A and
Subclade II-B is estimated within the last 1.25 mya followed by radiation within Subclade II-A and Subclade II-B was estimated at 0.65 and 0.85 mya, respectively.

| Diversity and haplotype analyses
Seven genetic indices were estimated from both the cytb and D loop dataset (only U.S. individuals; see Table 4). These included: TA B L E 2 (Continued) TA B L E 3 Average genetic distances of D-loop sequences estimated using the Tamura-Nei model of evolution (Tamura & Nei, 1993)

| Characterization of PRNP exon 3
DNA sequences from exon 3 of the PRNP gene were obtained from 10 individual aoudad revealed that all sequences were monomorphic.
Translation of nucleotides to amino acids revealed that aoudad possessed the signature genotype of A136, R154, and Q171 (Table 6), which is the most common genotype among sheep and goats (Goldmann, 2008).

| DISCUSS ION
It is important to note that data used in this research were obtained from three independent studies: cytb only (Derouiche et al., 2020), D loop only (Stipoljev et al., 2021), and cytb and D loop combined

TA B L E 3 (Continued)
F I G U R E 6 Time-calibrated phylogenetic tree modified from that depicted in Figure 42 with the superimposition of results from the BEAST analysis (version 2.6.1, Bouckaert et al., 2014)   Genetic divergence values between Clades I and II were 5.12% and 13.88% for cytb and D loop, respectively. Genetic distances obtained from the cytb gene (Table 2) indicated that the level of divergence between Clades I and II was much higher compared to values reported for other closely-related subspecies of bovids (x = 1.8%; e.g., Rupicapra rupicapra, Rupicapra pyrenaica, and Budorcas taxicolor, respectively). In addition, the genetic divergences of D loop observed between the two clades of aoudad (I and II) indicated a high level of genetic divergence compared to subspecies of Rupicapra (x = 8.04%; Table 3). The high levels of genetic divergence detected between Clades I and II indicates a magnitude of genetic divergence typically distinguishing subspecies of mammals (Baker & Bradley, 2006;Bradley & Baker, 2001).  (Tajima, 1989). In contrast, Fawcett WMA and Love Creek were identified by a positive Tajima's D, which may be indicative of balancing selection and population contraction (Tajima, 1989). For the D loop dataset, two populations (Fawcett WMA and Clade II as a whole) were identified by a test of Tajima's D as statistically significant ( Table 5). The negative values of Tajima's D for these two populations may be interpreted as population expansion as a result of a recent selective sweep (Tajima, 1989). Based on the available information garnered from translocation records, state agencies, and scientific publications (Barrett, 1980;Mungall & Sheffield, 1994;Ogren, 1959Ogren, , 1965Simpson & Krysl, 1981), it appears that the first translocated aoudad (most likely A.

Divergence dating analyses indicated that the
l. lervia, Cassinello, 1998;Gray, 1985;Ogren, 1965) in the United States initially were imported from European zoos to zoological parks and ultimately to private ranches. Sources indicated that the first free-ranging population in the United States was established on the Hearst Ranch (San Lucia Range, California) circa 1925 (Barrett, 1980) most likely sourced from the Fleishhacker Zoo (now the San Francisco Zoo; Mungall & Sheffield, 1994); unfortunately, importation records and other forms of documentation were unavailable and consequently were not useful in establishing the country of origin.
Escapees from the Hearst Ranch were thought to have established the contemporary population that currently is restricted to the San Lucia Range. Other individuals from the Hearst Ranch were used to establish zoo populations in California (San Francisco Zoo, Barrett, 1980;San Diego Zoo, Mungall & Sheffield, 1994 (Morrison, 1980) and Palo Duro Canyon in 1957and 1958(DeArment, 1971Mungall & Sheffield, 1994). In fact, the genetic data presented herein suggest that the most common haplogroup  (Simpson & Krysl, 1981) for the remaining two haplogroups restricted to Texas, the second-most common haplogroup (Clade I) appears to be representative of A. l. sahariensis or A. l. blainei. The difficulty in assigning subspecific origin (see Figure 3) to these samples stems from the fact that genetic data presented in Derouiche et al. (2020) seem to suggest that the Texas samples should be representative  Derouiche et al., 2020;Stipoljev et al., 2021). Further, Alados et al. (1988) and Castelló (2016) propose a wider distribution for A. l. blainei in the Ennedi and Uweinat mountains in northeast Chad, the native range of this subspecies currently is estimated to occur solely in the Red Hills of east Sudan (Cassinello, 2013;Nimir, 1997

| Concerns surrounding competition and extinction of sympatric haplogroups
Given the high level of genetic divergence (two major haplogroups -Clades I and II, Figures 2 and 3), it may be prudent to monitor populations to prevent homogenization of subspecies and loss of genetic and morphologic variation in light of their conservation status in their native range. First, it is noteworthy, that in Texas, the two major haplogroups (Clade I and Subclade II-B) were sympatric at six localities (Figure 3: 3,5,15,16,17,and 18). At these sites, haplogroup II-B was always present at a greater frequency, ranging from 57% to 91% (see Table 4). Further, two other localities ( On another note, aoudad exist in sympatry with native and charismatic Texas desert bighorn sheep in the Trans-Pecos Region and therefore the possibility for interspecies competition and disease transmission present concerns to the long-term management of these species (Barrett, 1967;Seegmiller & Simpson, 1979;Simpson & Krysl, 1981;Simpson et al., 1978). For example, throughout much of Texas and the desert southwest, water is a limited resource and results in the congregation of both species at natural and manmade water sources. Similarly, both species consume similar forage (Seegmiller & Simpson, 1979;Simpson et al., 1978) and may compete for this resource as well. Given the disparity in population size (>30,000 aoudad to 1500 bighorn sheep; F. Hernández, TPWD, personal communication), free-ranging populations of aoudad in the Trans-Pecos region will continue to increase as a result of year-round reproduction, high levels of recruitment, access to quality habitat and water resources, and minimal hunting pressure; thereby, outcompeting native bighorn sheep and possibly contribute to their extirpation (Barrett, 1967;Seegmiller & Simpson, 1979;Simpson & Krysl, 1981;Simpson et al., 1978). Further, the prion genotype, which confers average susceptibility to diseases such as scrapie (Goldmann, 2008), was detected in aoudad individuals examined in this study ( Table 6). Although the risk of prion transmission may be low, bighorn sheep and other native ungulates (e.g., mule deer, white-tailed deer, elk, and others) in this region may be at risk. For these and other reasons, population control of aoudad may become necessary to counter competition and disease transmission.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.